#include <algorithm>
#include <Geex/basics/assert.h>
#include "packer.h"


namespace Polygon_particle
{
	Packer* Packer::instance_ = 0 ;
	
	Packer::Packer(const std::string& polygon_fn, const std::string& debug_fn, double coverage_rate) : rand_gen(std::time(0))
	{
		gx_assert(instance_ == 0) ;
		instance_ = this ;
		std::srand(std::time(0));

		max_rate = coverage_rate;

		boost::random::mt19937 rand_gen(std::time(0));
		std::string ext = Geex::FileSystem::extension(polygon_fn);
		read_polygon_set(polygon_fn, std::back_inserter(polygon_set), ext);
		
		/////////////////////	Safe to remove	//////////////////////////
		// debug
		if (debug_fn.empty())
		{
			std::rotate(polygon_set.begin(), polygon_set.begin()+10, polygon_set.end());
	
			for (unsigned int i = 10; i < polygon_set.size(); i++)
				delete polygon_set[i];
			polygon_set.resize(10);
		}

		// debug file
		if (!debug_fn.empty())
		{
			std::ifstream db_file(debug_fn.c_str());
			while (db_file)
			{
				iter_trans.push_back(std::list<Parameter>());
				unsigned int i;
				for (i = 0; i < nb_polytopes(); i++)
				{
					double theta, x, y;
					if (! (db_file >> theta) ) break;
					if (! (db_file >> x) ) break;
					if (! (db_file >> y) ) break;
					iter_trans.back().push_back(Parameter(1.0, theta, x, y));
				}
				if (i < nb_polytopes())		iter_trans.pop_back();
			}
			mirror_polygons.resize(nb_polytopes());
			mirror_unions.resize(nb_polytopes());
			
			//////////////////////////////////////////////////////////////////
		}
		select_no = nb_polytopes();
		std::cout << "Specified coverage rate: "<<max_rate<<std::endl;

		double max_width = -std::numeric_limits<double>::max(), max_height = -std::numeric_limits<double>::max();
		total_volume = 0.0;
		for (unsigned int i = 0; i < polygon_set.size(); i++)
		{
			Packing_polygon::Iso_rectangle_2 bdbox = polygon_set[i]->bbox();
			max_width = std::max(bdbox.xmax() - bdbox.xmin(), max_width);
			max_height = std::max(bdbox.ymax() - bdbox.ymin(), max_height);
			if (polygon_set[i]->is_clockwise_oriented())
				polygon_set[i]->reverse_orientation();

			double vol = polygon_set[i]->volume();
			total_volume += vol;

		}
		max_width *= 3.0;
		rbox = new Rectangular_box(0.0, 0.0, max_width, total_volume/max_rate/max_width);
		cus.resize(nb_polytopes());

		///////////////////////	Initialization for energy method	//////////////////////////////////////
		if (debug_fn.empty())
		{
			max_height = -std::numeric_limits<double>::max();
			Bin_packing_2 bp(*rbox);
		
			for (unsigned int i = 0; i < polygon_set.size(); i++)
			{
				double x, y;
				if (bp.put_one_bin(polygon_set[i]->bbox(), x, y))
					polygon_set[i]->translate(x, y);
				max_height = std::max(polygon_set[i]->bbox().ymax(), max_height);
				//polygon_set[i]->update_skeleton();
				//cus[i] = polygon_set[i]->medial_circles();
			}
			
		}
		rbox->lower(rbox->y_max() - max_height);
		//////////////////////////////////////////////////////////////////////////////////////////////////
		compute_bounding_circles();
	
		std::cout<<nb_polytopes()<<" polygons are to be packed...\n";
	}
	Packer::~Packer() 
	{
		delete rbox;
		for (std::vector<Packing_polygon*>::iterator it = polygon_set.begin(); it != polygon_set.end(); ++it)
			delete *it;
		polygon_set.clear();
		instance_ = 0 ; 
	}

	// debug
	void Packer::trace()
	{
		static std::list< std::list<Parameter> >::iterator it = iter_trans.begin();
		if (it != iter_trans.end())
		{
			std::list<Parameter>::iterator one_trans = it->begin();
			for (unsigned int i = 0; i < nb_polytopes(); i++, ++one_trans)
			{
				mirror_polygons[i] = *polygon_set[i];
				mirror_polygons[i].transform(*one_trans);
				mirror_unions[i] = cus[i];
				mirror_unions[i].transform(one_trans->theta, one_trans->tx, one_trans->ty);
			}
			++it;
		}
		else 
			std::cout<<"End\n";
	}
	void Packer::hand_select(double mouse_x, double mouse_y)
	{
		Packing_polygon::Point mouse_pnt(mouse_x, mouse_y);
		for (unsigned int i = 0; i < nb_polytopes(); i++)
		{
			if (!polygon_set[i]->has_on_unbounded_side(mouse_pnt))
			{
				select_no = i;
				break;
			}
		}
	}
	void Packer::hand_move(double mouse_x, double mouse_y)
	{
		if (select_no >= nb_polytopes())
			return;
		Packing_polygon::Point mouse_pnt(mouse_x, mouse_y);
		Packing_polygon::Point cent = polygon_set[select_no]->centroid();
		Packing_polygon::Vector v(cent, mouse_pnt);
		polygon_set[select_no]->translate(v.x(), v.y());
		cus[select_no].transform(0.0, v.x(), v.y());
		select_no = nb_polytopes();
	}
	void Packer::print_polygons()
	{
		std::ofstream ofs("init_polygons.poly");
		ofs << nb_polytopes() << std::endl;
		for (unsigned int i = 0; i < nb_polytopes(); i++)
		{
			ofs << polygon_set[i]->size() << std::endl;
			for (unsigned int j = 0; j < polygon_set[i]->size(); j++)
			{
				ofs << polygon_set[i]->vertex(j).x() << ' ';
				ofs << polygon_set[i]->vertex(j).y() << '\n';
			}
		}
	}
	void Packer::compute_bounding_circles()
	{
		const unsigned int N = 20; // referential circle number
		cus.resize(nb_polytopes());
		for (unsigned int i = 0; i < polygon_set.size(); i++)
		{
			// decide the initial offset
			double area = std::fabs(polygon_set[i]->area());
			double offset = 0.2*area/(N*0.5*boost::math::constants::pi<double>());
			offset = std::sqrt(offset);
			std::vector< boost::shared_ptr<CGAL_polygon> > offset_polygon = CGAL::create_interior_skeleton_and_offset_polygons_2(offset, *dynamic_cast<CGAL_polygon*>(polygon_set[i]));
			while (offset_polygon.size() != 1)
			{
				offset *= 0.98;
				offset_polygon = CGAL::create_interior_skeleton_and_offset_polygons_2(offset, *dynamic_cast<CGAL_polygon*>(polygon_set[i]));
			}
			// sampling
			//unsigned int pre = offset_polygon.front()->size() - 1;
			//unsigned int sz = offset_polygon.front()->size();
			// store the circles at the vertices first
			for (unsigned int j = 0; j < offset_polygon.front()->size(); j++)
			{
				//Packing_polygon::Vector v0 = offset_polygon.front()->edge(j).to_vector(), v1 = - offset_polygon.front()->edge(pre).to_vector();
				//double cos_2theta = v0*v1/CGAL::sqrt(v0.squared_length()*v1.squared_length());
				//cos_2theta = std::max(cos_2theta, -1.0);
				//cos_2theta = std::min(cos_2theta, 1.0);
				//double theta = 0.5*std::acos(cos_2theta);
				//double sin_theta = std::sin(theta);
				//double r = offset/sin_theta;
				double r = offset;
				Packing_polygon::Point src = offset_polygon.front()->edge(j).source()/*, tgt = offset_polygon.front()->edge(j).target()*/;
				cus[i].push_back(Particle_circle(src, r*r));
			}

			for (unsigned int j = 0; j < offset_polygon.front()->size(); j++)
			{
				Packing_polygon::Segment_2 edge = offset_polygon.front()->edge(j);
				double edge_len = CGAL::sqrt(edge.squared_length());
				//int n = (edge_len - cus[i].circle_at(j).radius() - cus[i].circle_at((j+1)%sz).radius())/(2*offset);
				int n = (edge_len - 2*offset)/(2*offset);
				if (n <= 0)	continue;
				//unsigned int n = edge_len/(2*offset);
				n = std::max<unsigned int>(n, 1);
				double step = edge_len/(n+1);
				Packing_polygon::Vector v = edge.to_vector()/edge_len;
				Packing_polygon::Point c = offset_polygon.front()->edge(j).source() + v*step;
				for (unsigned int k = 0; k < n; k++, c = c + v*step)	
				{
					Particle_circle cir(c, offset*offset);
					// ASSUME all the polygons are in the counterclockwise order
					Particle_circle::Vector_2 nm = edge.to_vector().perpendicular(CGAL::COUNTERCLOCKWISE);
					nm = nm/CGAL::sqrt(nm.squared_length());
					cir.orientation_vector(nm);
					cus[i].push_back(cir);					
				}
			}
#if 0
			unsigned int pre = polygon_set[i]->size()-1;
			unsigned int sz = polygon_set[i]->size();
			for (unsigned int j = 0; j < sz; j++, pre = (pre+1)%sz)
			{
				Packing_polygon::Vector v0 = polygon_set[i]->edge(j).to_vector(), v1 = -polygon_set[i]->edge(pre).to_vector();
				v0 = v0/CGAL::sqrt(v0.squared_length());
				v1 = v1/CGAL::sqrt(v1.squared_length());
				double cos_2theta = v0*v1;
				cos_2theta = std::max(cos_2theta, -1.0);
				cos_2theta = std::min(cos_2theta, 1.0);
				double theta = 0.5*std::acos(cos_2theta);
				unsigned int n = boost::math::constants::pi<double>()*0.5/theta+1;
				double sin_theta = std::sin(theta), cos_theta = std::cos(theta);
				Packing_polygon::Vector bisec = v0 + v1;
				bisec = bisec/CGAL::sqrt(bisec.squared_length());
				double arc_len = offset/sin_theta - offset;
				double step = arc_len/n;
				Packing_polygon::Segment_2 s = polygon_set[i]->edge(j);
				for (unsigned int k = 0; k < n; k++)
				{
					Packing_polygon::Point c = s.source() + k*step*bisec;
					double sr = CGAL::squared_distance(c, s.supporting_line());
					if (sr > 0.0)
						cus[i].push_back(Particle_circle(c, sr));
				}
			}
#endif
		}
	}

	void Packer::pack(unsigned int nb_trials)
	{

		//for (unsigned int i = 0; i < nb_trials; i++)
		//{
		//	double e = energy_minimize();
		//	if (e < min_energy)
		//	{
		//		min_energy = e;
		//		best_config.assign(circles.begin(), circles.end());
		//	}
		//	glut_viewer_redraw();
		//}
		//std::cout << "Minimal energy found = "<<min_energy<<std::endl;
		//circles.assign(best_config.begin(), best_config.end());
		//glut_viewer_redraw();
	}
	
	double Packer::energy_minimize()
	{
		static double opti_time = 0.0;
		static unsigned int opti_nb = 0;

//		CGAL::Timer timer;
		System_potential_energy *spe = new Circle_union_system(cus, *dynamic_cast<Rectangular_box*>(rbox));
		double *trans = new double[cus.size()*3];
		std::fill(trans, trans + cus.size()*3, 0.0);

		//timer.start();

		int res = spe->optimize(trans, false);
		double e = std::numeric_limits<double>::max();	

		if (res >= -400)
		{
			for (unsigned int i = 0; i < cus.size(); i++)
			{
				//std::cout<<'['<<trans[i*3]<<','<<trans[3*i+1]<<','<<trans[3*i+2]<<"], ";
				cus[i].transform(trans[i*3], trans[3*i+1], trans[3*i+2]);	
				polygon_set[i]->transform(Parameter(1.0, trans[i*3], trans[i*3+1], trans[i*3+2]));
			}

			//std::cout<<"Try to eliminate overlap...\n";
			//std::fill(trans, trans + cus.size()*3, 0.0);
			//res = spe->optimize(trans, true);
			//if (res >= -400)
			//{
			//	// after query, we do update
			//	double ge = spe->gravity_energy_at(trans);
			//	e = spe->system_energy_at(trans);
			//	std::cout<<"Potential energy = "<<e<<", gravity energy = "<<ge<<", LJ energy = "<<e-ge<<std::endl;
			//	for (unsigned int i = 0; i < cus.size(); i++)
			//	{
			//		cus[i].transform(trans[i*3], trans[3*i+1], trans[3*i+2]);	
			//		polygon_set[i]->transform(Parameter(1.0, trans[i*3], trans[i*3+1], trans[i*3+2]));
			//	}
			//}
			//else
			//{
			//	std::cout<<"Failed to eliminate overlap, with error code = "<<res<<std::endl;
			//}
		}
		else
		{
			std::cout<<"Minimization failed, error code = "<<res<<std::endl;
		}

		//timer.stop();
		//opti_time += timer.time();
		//opti_nb++;
		//std::cout<<"Average optimization time = "<<opti_time/opti_nb<<std::endl;
		delete[] trans;
		delete spe;
		return e;
	}
	void Packer::perturb()
	{
		//boost::random::uniform_real_distribution<double> uni_distr(0.0, boost::math::constants::two_pi<double>());
		//boost::variate_generator< boost::random::mt19937&, boost::random::uniform_real_distribution<double> > gen(rand_gen, uni_distr);
		//for (unsigned int i = 0 ; i < circles.size(); i++)
		//{
		//	double rot_angle = gen();
		//	Particle_circle::Vector_2 dir(std::cos(rot_angle), std::sin(rot_angle));
		//	double len = (gen()/10)*circles[i].radius();
		//	circles[i].translate(len*dir.x(), len*dir.y());
		//}
	}
	void Packer::shuffle()
	{
	}


	void Packer::get_bbox(Geex::real& x_min, Geex::real& y_min, Geex::real& z_min, Geex::real& x_max, Geex::real& y_max, Geex::real& z_max)
	{
		x_min = rbox->x_min(); x_max = rbox->x_max();
		y_min = rbox->y_min(); y_max = rbox->y_max();
		z_min = rbox->z_min(); z_max = rbox->z_max();
	}
}
